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^ Abstract 

^H High dimensional time series are endemic in applications of machine learning such as robotics (sensor 

^H data), computational biology (gene expression data), vision (video sequences) and graphics (motion capture 

^— > data). Practical nonlinear probabilistic approaches to this data are required. In this paper we introduce the 

1^ variational Gaussian process dynamical system. Our work builds on recent variational approximations for 

C^ Gaussian process latent variable models to allow for nonlinear dimensionality reduction simultaneously with 

learning a dynamical prior in the latent space. The approach also allows for the appropriate dimensionality of 

I the latent space to be automatically determined. We demonstrate the model on a human motion capture data 

^ set and a series of high resolution video sequences. 

in 
oo 
a^ 1 Introduction 

f- — . Nonlinear probabilistic modeling of high dimensional time series data is a key challenge for the machine learn- 

^^ ing community. A standard approach is to simultaneously apply a nonlinear dimensionality reduction to the 

^■^ data whilst governing the latent space with a nonlinear temporal prior. The key difficulty for such approaches is 

^7; that analytic marginalization of the latent space is typically intractable. Markov chain Monte Carlo approaches 

^ can also be problematic as latent trajectories are strongly correlated making efficient sampling a challenge. One 

promising approach to these time series has been to extend the Gaussian process latent variable model [1,2] 
with a dynamical prior for the latent space and seek a maximum a posteriori (MAP) solution for the latent 
C^ points [3, 4, 5]. Ko and Fox [6] further extend these models for fully Bayesian filtering in a robotics setting. We 

refer to this class of dynamical models based on the GP-LVM as Gaussian process dynamical systems (GPDS). 
However, the use of a MAP approximation for training these models presents key problems. Firstly, since the 
latent variables are not marginalised, the parameters of the dynamical prior cannot be optimized without the 
risk of overfitting. Further, the dimensionality of the latent space cannot be determined by the model: adding 
further dimensions always increases the likelihood of the data. In this paper we build on recent developments 
in variational approximations for Gaussian processes [7, 8] to introduce a variational Gaussian process dynami- 
cal system (VGPDS) where latent variables are approximately marginalized through optimization of a rigorous 
lower bound on the marginal likelihood. As well as providing a principled approach to handling uncertainty in 
the latent space, this allows both the parameters of the latent dynamical process and the dimensionality of the 
latent space to be determined. The approximation enables the application of our model to time series containing 
millions of dimensions and thousands of time points. We illustrate this by modeling human motion capture data 
and high dimensional video sequences. 



^ 



2 The Model 

Assume a multivariate times series dataset {yn,tn}^=i, where y^ G M^ is a data vector observed at time 
tn G R+. We are especially interested in cases where each y^ is a high dimensional vector and, therefore, 
we assume that there exists a low dimensional manifold that governs the generation of the data. Specifically, a 
temporal latent function x(t) G MP (with Q <C I^), governs an intermediate hidden layer when generating the 
data, and the di\\ feature from the data vector y^ is then produced from x^ = x(tn) according to 

2/nd =/d(Xn) +end , e^d -- A/'(0, /3"^), (1) 

where fd (x) is a latent mapping from the low dimensional space to di\\ dimension of the observation space 
and P is the inverse variance of the white Gaussian noise. We do not want to make strong assumptions about 
the functional form of the latent functions (x, f).^ Instead we would like to infer them in a fully Bayesian 
non-parametric fashion using Gaussian processes [9] . Therefore, we assume that x is a multivariate Gaussian 
process indexed by time t and f is a different multivariate Gaussian process indexed by x, and we write 

Xq(t)r^gV({),k,{U,tj)), g = l,...,Q, (2) 

/rf(x)^anO,fe/(x„x,)), d=l,...,D. (3) 

Here, the individual components of the latent function x are taken to be independent sample paths drawn from 
a Gaussian process with covariance function k^ {ti ^tj). Similarly, the components of f are independent draws 
from a Gaussian process with covariance function kf{^i^^j). These covariance functions, parametrized by 
parameters Ox and Of respectively, play very distinct roles in the model. More precisely, kx determines the 
properties of each temporal latent function Xq{t). For instance, the use of an Ornstein-Uhlbeck covariance 
function yields a Gauss-Markov process for Xq(t), while the squared-exponential kernel gives rise to very 
smooth and non-Markovian processes. In our experiments, we will focus on the squared exponential covariance 
function (RBF), the Matem 3/2 which is only once differentiable, and a periodic covariance function [9, 10] 
which can be used when data exhibit strong periodicity. These kernel functions take the form: 



^x(rbf) {U-, tj) = CTj-bfe (^^?) , kx(mat) {U,tj) = CT, 



2 



V^\U-tj\\ z^^l^izhl 







^x(per) {U,tj) =(T^^,e 2 h . (4) 



The covariance function kf determines the properties of the latent mapping f that maps each low dimensional 
variable x^ to the observed vector y^. We wish this mapping to be a non-linear but smooth, and thus a suitable 
choice is the squared exponential covariance function 

kf (xi,x,-) = al^e-'^ T.U-'^{-'^^-=.)\ (5) 

which assumes a different scale Wq for each latent dimension. This, as in the variational Bayesian formulation 
of the GP-LVM [8], enables an automatic relevance determination procedure (ARD), i.e. it allows Bayesian 
training to "switch off" unnecessary dimensions by driving the values of the corresponding scales to zero. 

The matrix Y G R^^^ will collectively denote all observed data so that its nth row corresponds to the 
data point y^. Similarly, the matrix F G R^^^ will denote the mapping latent variables, i.e. fnd — fdi'^n)^ 
associated with observations Y from (1). Ana usly, X G R^^^ will store all low dimensional latent variables 
Xnq = Xq{tn)- Further, we will refer to columns of these matrices by the vectors y^^, f^^, x^ G R^. Given the 
latent variables we assume independence over the data features, and given time we assume independence over 
latent dimensions to give 

D Q 

p{Y,F,X\t) =p{Y\F)p{F\X)p{X\t) = l[p{yd\U)p{U\X) l[p{^q\t)^ (6) 

d=l q=l 



^To simplify our notation, we often write x instead of x(t) and f instead of f (x). Later we also use a similar convention for the kernel 
functions by often writing them as kf and kx . 



where t G M^ and p{yd\^d) is a Gaussian likelihood function term defined from (1). Further, p(f(^|X) is a 
marginal GP prior such that 

p{U\X)=M{U\O^Knn). (7) 

where K^n = kf{X^ X) is the covariance matrix defined by the kernel function kf and similarly p(xg|t) is the 
marginal GP prior associated with the temporal function Xq{t), 

p(x,|t)=A/'(x,|0,ir,), (8) 

where Kf = /Ca^(t, t) is the covariance matrix obtained by evaluating the kernel function k^ on the observed 
times t. 

Bayesian inference using the above model poses a huge computational challenge as, for instance, marginal- 
ization of the variables X, that appear non-linearly inside the kernel matrix K^n^ is troublesome. Practical 
approaches that have been considered until now (e.g. [5, 3]) marginalise out only F and seek a MAP solution 
for X. In the next section we describe how efficient variational approximations can be applied to marginalize 
X by extending the framework of [8]. 

2.1 Variational Bayesian training 

The key difficulty with the Bayesian approach is propagating the prior density p{X\t) through the nonlinear 
mapping. This mapping gives the expressive power to the model, but simultaneously renders the associated 
marginal likelihood, 

p(Y\t) = I p{Y\F)p{F\X)p{X\t)dXdF, (9) 

intractable. We now invoke the variational Bayesian methodology to approximate the integral. Following a 
standard procedure [1 1], we introduce a variational distribution q{Q) and compute the Jensen's lower bound Ty 
on the logarithm of (9), 

^„(,,.) = /,(e,iog'ia£MQS£(a),xdf, (10, 

where denotes the model's parameters. However, the above form of the lower bound is problematic becauce 
X (in the GP term p{F\X)) appears non-linearly inside the kernel matrix Knn making the integration over X 
difficult. As shown in [8], this intractability is removed by applying the "data augmentation" principle. More 
precisely, we augment the joint probability model in (6) by including M extra samples of the GP latent mapping 
f , known as inducing points, so that u^ G M^ is such a sample. The inducing points are evaluated at a set of 
pseudo-inputs X G R^^^. The augmented joint probability density takes the form 

D 

p(r,F,/7,x,x|t) = np(yd|fdMfd|u,,x)p(u,|x)p(x|t), (ii) 

where p{ud\X) is 3. zero-mean Gaussian with a covariance matrix Kmm constructed using the same function 
as for the GP prior (7). By dropping X from our expressions, we write the augmented GP prior analytically 
(see [9]) as 

p(frf|urf, X) = A/" (frfli^ATMi^MMUrf, Knn - KnmKI^^Kmn) ' (12) 

A key result in [8] is that a tractable lower bound (computed analogously to (10)) can be obtained through the 
variational density 

D 

q{0) = q{F, U,X) = q{F\U,X)q{U)q{X) = [] p(f^|u^, X)^(u^)^(X), (13) 

where g(X) = Hg^i -^ (^glMg^ ^q) ^^^ Qi'^d) is an arbitrary variational distribution. Titsias and Lawrence [8] 
assume full independence for q{X) and the variational covariances are diagonal matrices. Here, in contrast, the 



posterior over the latent variables will have strong correlations, so Sq is taken to be a A^ x A^ full co variance 
matrix. Optimization of the variational lower bound provides an approximation to the true posterior p{X\Y) 
by q{X). In the augmented probability model, the "difficult" term p{F\X) appearing in (10) is now replaced 
with (12) and, eventually, it cancels out with the first factor of the variational distribution (13) so that F can be 
marginalised out analytically. Given the above and after breaking the logarithm in (10), we obtain the final form 
of the lower bound (see supplementary material for more details) 

J^,{q, e)=P,- KL(g(X) II p(X|t)), (14) 

with Py = J q{X) \ogp{Y\F)p{F\X) dXdF. Both terms in (14) are now tractable. Note that the first of 
the above terms involves the data while the second one only involves the prior. All the information regarding 
data point correlations is captured in the KL term and the connection with the observations comes through the 
variational distribution. Therefore, the first term in (14) has the same analytical solution as the one derived in 
[8]. (14) can be maximized by using gradient-based methods^. However, when not factorizing q{X) across data 
points yields 0{N'^) variational parameters to optimize. This issue is addressed in the next section. 

2.2 Reparametrization and Optimization 

The optimization involves the model parameters 6 = (/3, ^/, 6t), the variational parameters {)Ug, Sq}^^^ from 
q{X) and the inducing points^ X. 

Optimization of the variational parameters appears challenging, due to their large number and the correla- 
tions between them. However, by reparametrizing our O {N'^) variational parameters according to the frame- 
work described in [12] we can obtain a set of 0(7V) less correlated variational parameters. Specifically, we first 
take the derivatives of the variational bound (14) w.r.t. Sq and /ig and set them to zero, to find the stationary 
points, 

Sq = {K-^ + Aq)'^ and flq = Ktflq, (15) 

where A^ = —2 — ^y' ■' is sl N x N diagonal, positive matrix and flq = |^ is a A/"— dimensional vector. 
The above stationary conditions tell us that, since Sq depends on a diagonal matrix A^, we can reparametrize it 
using only the X— dimensional diagonal of that matrix, denoted by Xq. Then, we can optimise the 2{Q x N) 
parameters {Xq, flq) and obtain the original parameters using (15). 

2.3 Learning from Multiple Sequences 

Our objective is to model multivariate time series. A given data set may consist of a group of independent ob- 
served sequences, each with a different length (e.g. in human motion capture data several walks from a subject). 
Let, for example, the dataset be a group of S independent sequences (Y^^"^ , ..., Y^^"^). We would like our model 
to capture the underlying commonality of these data. We handle this by allowing a different temporal latent 
function for each of the independent sequences, so that X'^^^ is the set of latent variables corresponding to the 
sequence s. These sets are a priori assumed to be independent since they correspond to separate sequences, 
i.e. p (A^^\ A^^\ ..., A^*^^) = ris^i P(^^^^)' where we dropped the conditioning on time for simplicity. This 
factorisation leads to a block-diagonal structure for the time covariance matrix Kt, where each block corre- 
sponds to one sequenece. In this setting, each block of observations Y^^"^ is generated from its corresponding 
A*^*^ according to Y'^^^ = F*^*^ + e, where the latent function which governs this mapping is shared across all 
sequences and e is Gaussian noise. 



^See supplementary material for more detailed derivation of (14) and for the equations for the gradients. 

^We will use the term "variational parameters" to refer only to the parameters of q{X) although the inducing points are also variational 
parameters. 



3 Predictions 

Our algorithm models the temporal evolution of a dynamical system. It should be capable of generating com- 
pletely new sequences or reconstructing missing observations from partially observed data. For generating novel 
sequence given training data the model requires a time vector t* as input and computes a density p{Y^\Y^ t, t*). 
For reconstruction of partially observed data the time stamp information is additionally accompanied by a par- 
tially observed sequence Y^ G M^*x^p from the whole Y^ = (Kf , Y^), where p and m are set of indices 
indicating the present (i.e. observed) and missing dimensions of Y^ respectively, so that pU m = {1, . . . , D}. 
We reconstruct the missing dimensions by computing the Bayesian predictive distribution p{Y^\Yi^^ F, t^, t). 
The predictive densities can also be used as estimators for tasks like generative Bayesian classification. Whilst 
time stamp information is always provided, in the next section we drop its dependence to avoid notational 
clutter. 

3.1 Predictions Given Only the Test Time Points 

To approximate the predictive density, we will need to introduce the underlying latent function values F^ G 
R^* ^^ (the noisy-free version of Y^) and the latent variables X^ G R^* x^. We write the predictive density as 

p{Y,\Y) = fp{Y,,F,,X,\Y,,Y)dF,dX, = f p{Y,\F,)p{F,\X,,Y)p{X,\Y)dF,dX,. (16) 

The term p{F^\X^^Y) is approximated by the variational distribution 

q{F,\X,) = / n P(f*,d|ud,X*)g(urf)durf = II g(f*,d|X*), (17) 

-^ deD deD 

where q{U,d\X^) isa. Gaussian that can be computed analytically, since in our variational framework the optimal 
setting for q{ud) is also found to be a Gaussian (see suppl. material for complete forms). As for the term 
p{X^\Y) in eq. (16), it is approximated by a Gaussian variational distribution q{X^), 

Q Q . Q 

Q{^*) = n ^(^*'^) = n / P(x*,g|xJ(7(Xg)dXg = Yl {P{^*,q\^q))qi^^) ^ (1^) 

q=l q=l ^ q=l 

where p(x^ q|xg) is a Gaussian found from the conditional GP prior (see [9]) and q{X) is also Gaussian. We 
can, thus, work out analytically the mean and variance for (18), which turn out to be: 

/ia,^^ = K^Npq (19) 

var(x*,J = K^^ - K^N{Kt + A^^)"^i^Ar*. (20) 

where K^n = ^a?(t*,t), K^n = Kjj^ and K^^ = kx{t^^t^). Notice that these equations have exactly the 
same form as found in standard GP regression problems. Once we have analytic forms for the posteriors in (16), 
the predictive density is approximated as 

p{Y,\Y) = f p{Y,\F,)q{F,\X,)q{X,)dF,dX, = f p{Y,\F,) {q{F,\X,))^^^^^dF,, (21) 

which is a non-Gaussian integral that cannot be computed analytically. However, following the same argument 
as in [9, 13], we can calculate analytically its mean and covariance: 

E(F*) = B^^l (22) 

Cov(F*) = B'^ (^* - ^*(^t)^) B^%I- Tr [(i^^^^ - {Kmm + /3^2)"') ^2] ^, (23) 

where B = /3 (i^MM + /3^2)"' ^7^, ^ = (%(X*,X*)), ^* = (Km.) and ^J = {Km.K.m}- All 
expectations are taken w.r.t. q{X^) and can be calculated analytically, while Km^ denotes the cross-co variance 
matrix between the training inducing inputs X and X^. The ^ quantities are calculated analytically (see suppl. 
material). Finally, since Y^ is just a noisy version of F^, the mean and covariance of (21) is just computed as: 

E(i;) =E(F*)andCov(i;) =Cov(F*)+/3-^/Ar,. 



3.2 Predictions Given the Test Time Points and Partially Observed Outputs 

The expression for the predictive density p{Y^\Yi^^ Y) is similar to (16), 

p(rr ii?, Y) = j p{Yr\Fr)p{Fr\x,, i?, FMx.iyr, r)dFrdx„ (24) 

and is analytically intractable. To obtain an approximation, we firstly need to apply variational inference and 
approximate p{X^\Yi ^Y) with a Gaussian distribution. This requires the optimisation of a new variational 
lower bound that accounts for the contribution of the partially observed data Y^ . This lower bound approximates 
the true marginal likelihood p{Yi ^Y) and has exactly analogous form with the lower bound computed only on 
the training data Y . Moreover, the variational optimisation requires the definition of the variational distribution 
q{X^^X) which needs to be optimised and is fully correlated across X and X^. After the optimisation, the 
approximation to the true posterior p(X*|y?, Y) is given from the marginal q{X^). A much faster but less 
accurate method would be to decouple the test from the training latent variables by imposing the factorisation 
q{X^^X) = q{X)q{X^). This is not used, however, in our current implementation. 

4 Handling Very High Dimensional Datasets 

Our variational framework avoids the typical cubic complexity of Gaussian processes allowing relatively large 
training sets (thousands of time points, N). Further, the model scales only linearly with the number of dimen- 
sions D. Specifically, the number of dimensions only matters when performing calculations involving the data 
matrix Y . In the final form of the lower bound (and consequently in all of the derived quantities, such as gra- 
dients) this matrix only appears in the form YY^ which can be precomputed. This means that, when N ^ D, 
we can calculate YY^ only once and then substitute Y with the SVD (or Cholesky decomposition) of YY^ . In 
this way, we can work with an A^ x TV instead of an A" x D matrix. Practically speaking, this allows us to work 
with data sets involving millions of features. In our experiments we model directly the pixels of HD quality 
video, exploiting this trick. 

5 Experiments 

We consider two different types of high dimensional time series, a human motion capture data set consist- 
ing of different walks and high resolution video sequences. The experiments are intended to explore the 
various properties of the model and to evaluate its performance in different tasks (prediction, reconstruction, 
generation of data). Matlab source code for repeating the following experiments is available on-line from 

http : //staff www. dcs . shef . ac . uk/people/N. Lawrence /vargplvm/. 

5.1 Human Motion Capture Data 

We followed [14, 15] in considering motion capture data of walks and runs taken from subject 35 in the CMU 
motion capture database. We treated each motion as an independent sequence. The data set was constructed and 
preprocessed as described in [15]. This results in 2,613 separate 59-dimensional frames split into 31 training 
sequences with an average length of 84 frames each. 

The model is jointly trained, as explained in section 2.3, on both walks and runs, i.e. the algorithm learns 
a common latent space for these motions. At test time we investigate the ability of the model to reconstruct 
test data from a previously unseen sequence given partial information for the test targets. This is tested once 
by providing only the dimensions which correspond to the body of the subject and once by providing those that 
correspond to the legs. We compare with results in [15], which used MAP approximations for the dynamical 
models, and against nearest neighbour. We can also indirectly compare with the binary latent variable model 
(BLV) of [14] which used a slightly different data preprocessing. We assess the performance using the cumu- 
lative error per joint in the scaled space defined in [14] and by the root mean square error in the angle space 
suggested by [15]. Our model was initialized with nine latent dimensions. We performed two runs, once using 



the Matem covariance function for the dynamical prior and once using the RBF. From table 1 we see that the 
variational Gaussian process dynamical system considerably outperforms the other approaches. The appropriate 
latent space dimensionality for the data was automatically inferred by our models. The model which employed 
an RBF covariance to govern the dynamics retained four dimensions, whereas the model that used the Matem 
kept only three. The other latent dimensions were completely switched off by the ARD parameters. The best 
performance for the legs and the body reconstruction was achieved by the VGPDS model that used the Matem 
and the RBF covariance function respectively. 

Table 1 : Errors obtained for the motion capture dataset considering nearest neighbour in the angle space (NN) and in the 
scaled space(NN sc), GPLVM, BLV and VGPDS. CL / CB are the leg and body datasets as preprocessed in [14], L and B 
the corresponding datasets from [15]. SC corresponds to the error in the scaled space, as in Taylor et al. while RA is the 
error in the angle space. The best error per column is in bold. 



Data 


CL 


CB 


L 


L 


B 


B 


Error Type 


SC 


SC 


SC 


RA 


SC 


RA 


BLV 


n.7 


8.8 


- 


- 


- 


- 


NNsc. 


22.2 


20.5 


- 


- 


- 


- 


GPLVM (Q = 3) 


- 


- 


n.4 


3.40 


16.9 


2.49 


GPLVM (Q = 4) 


- 


- 


9.7 


3.38 


20.7 


2.72 


GPLVM (Q = 5) 


- 


- 


13.4 


4.25 


23.4 


2.78 


NNsc. 


- 


- 


13.5 


4.44 


20.8 


2.62 


NN 


- 


- 


14.0 


4.11 


30.9 


3.20 


VGPDS (RBF) 


- 


- 


8.19 


3.57 


10.73 


1.90 


VGPDS (Matem 3/2) 


- 


- 


6.99 


2.88 


14.22 


2.23 



5.2 Modeling Raw High Dimensional Video Sequences 

For our second set of experiments we considered video sequences. Such sequences are typically preprocessed 
before modeling to extract informative features and reduce the dimensionality of the problem. Here we work 
directly with the raw pixel values to demonstrate the ability of the VGPDS to model data with a vast number of 
features. This also allows us to directly sample video from the learned model. 

Firstly, we used the model to reconstruct partially observed frames from test video sequences'^. For the first 
video discussed here we gave as partial information approximately 50% of the pixels while for the other two we 
gave approximately 40% of the pixels on each frame. The mean squared error per pixel was measured to com- 
pare with the /c— nearest neighbour (NN) method, fork G (1, .., 5) (we only present the error achieved for the 
best choice of k in each case). The datasets considered are the following: firstly, the 'Missa' dataset, a standard 
benchmark used in image processing. This is 103,680-dimensional video, showing a woman talking for 150 
frames. The data is challenging as there are translations in the pixel space. We also considered an HD video of 
dimensionality 9 x 10^ that shows an artificially created scene of ocean waves as well as a 230, 400— dimensional 
video showing a dog running for 60 frames. The later is approximately periodic in nature, containing several 
paces from the dog. For the first two videos we used the Matern and RBF kernel respectively to model the dy- 
namics and interpolated to reconstruct blocks of frames chosen from the whole sequence. For the 'dog' dataset 
we constructed a compound kernel kx — /^^(rbf) + ^x (periodic) ^ where the RBF term is employed to capture any 
divergence from the approximately periodic pattern. We then used our model to reconstruct the last 7 frames ex- 
trapolating beyond the original video. As can be seen in table 2, our method outperformed NN in all cases. The 
results are also demonstrated visually in figure 1 and the reconstructed videos are available in the supplementary 
material. 



"^ 'Missa' dataset: cipr.rpi.edu. 'Ocean': cogfilms.com. 'Dog': fitfurlife.com. See details in supplementary. The logo appearing in the 
'dog' images in the experiments that follow, has been added with post-processing. 



Table 2: The mean squared error per pixel for VGPDS and NN for the three datasets (measured only in the missing inputs). 
The number of latent dimensions selected by our model is in parenthesis. 





Missa 




Ocean 




Dog 




VGPDS 


2.52 (Q = 


= 12) 


9.36 (Q = 


= 9) 


4.01 (Q = 


= 6) 


NN 


2.63 




9.53 




4.15 






(e) 



(f) 



(g) 



(h) 



Figure 1 : (a) and (c) demonstrate the reconstruction achieved by VGPDS and NN respectively for the most challenging 
frame (b) of the 'missa' video, i.e. when translation occurs, (d) shows another example of the reconstruction achieved by 
VGPDS given the partially observed image, (e) (VGPDS) and (f) (NN) depict the reconstruction achieved for a frame of 
the 'ocean' dataset. Finally, we demonstrate the ability of the model to automatically select the latent dimensionality by 
showing the initial lengthscales (fig: (g)) of the ARD kernel and the values obtained after training (fig: (h)) on the 'dog' data 
set. 



As can be seen in figure 1, VGPDS predicts pixels which are smoothly connected with the observed image, 
whereas the NN method cannot fit the predicted pixels in the overall context. 

As a second task, we used our generative model to create new samples and generate a new video sequence. 
This is most effective for the 'dog' video as the training examples were approximately periodic in nature. The 
model was trained on 60 frames (time-stamps [ti, teo]) and we generated the new frames which correspond to 
the next 40 time points in the future. The only input given for this generation of future frames was the time 
stamp vector, [tei, tioo]. The results show a smooth transition from training to test and amongst the test video 
frames. The resulting video of the dog continuing to run is sharp and high quality. This experiment demonstrates 
the ability of the model to reconstruct massively high dimensional images without blurring. Frames from the 
result are shown in figure 2. The full video is available in the supplementary material. 



6 Discussion and Future Work 



We have introduced a fully Bayesian approach for modeling dynamical systems through probabilistic nonlinear 
dimensionality reduction. Marginalizing the latent space and reconstructing data using Gaussian processes 
results in a very generic model for capturing complex, non-linear correlations even in very high dimensional 



Figure 2: The last frame of the training video (a) is smoothly followed by the first frame (b) of the generated video. A 
subsequent generated frame can be seen in (c). 

data, without having to perform any data preprocessing or exhaustive search for defining the model's structure 
and parameters. 

Our method's effectiveness has been demonstrated in two tasks; firstly, in modeling human motion capture 
data and, secondly, in reconstructing and generating raw, very high dimensional video sequences. A promising 
future direction to follow would be to enhance our formulation with domain- specific knowledge encoded, for 
example, in more sophisticated covariance functions or in the way that data are being preprocessed. Thus, we 
can obtain application-oriented methods to be used for tasks in areas such as robotics, computer vision and 
finance. 
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Appendix 
A Derivation of the variational bound 

We wish to approximate the marginal HkeHhood: 

p{Y\t) = Jp{Y,F,X\t)dXdF, (25) 

by computing a lower bound: 

Ml, 0)= Jq{ O) log ?^^2l^AXdF, (26) 

This can be achieved by first augmenting the joint probability density of our model with inducing inputs X 
along with their corresponding function values U: 

D 

p{Y,F,U,X,X\t) = \{p{ya\id)p(.U^d,X)p{ud\X)p{X\t) (27) 

where p{\id\X) = fld^i -^ {^d\^-, Kmm) • For simplicity, X is dropped from our expressions for the rest of 
this supplementary material. Note that after including the inducing points, p(fc^|uc^,X) remains analytically 
tractable and it turns out to be [9]): 

v{U\^d,X) = M (fd\KNMK^j^Ud,KNN - KnmK^^Kmn) • (28) 

We are now able to define a variational distribution g(B) which factorises as: For tractability we now define a 
variational density, q{Q)\ 

D 

q{e) = q{F,U,X) = q{F\U,X)q{U)q{X) = ^^(f.lu,, X)g(u,)g(X), (29) 

d=l 

where q{X) = YYq^iJ^ {^qlfJ^q^ Sq). Now, we return to (26) and replace the joint distribution with its aug- 
mented version (27) and the variational distribution with its factorised version (29): 

= / n rtfju.. A-)<,(u,)<,(A^) log n?^- 'Jff'lf^f:^'^''']^''™ dXdf 

p{id\ud, X)q{ud)q{X) log -^^ dXdF, 

d=i lld=i q{ud)q{X)) 

= T,-KL{q\\p), (30) 

with Pv = j q{X) \ogp{Y\F)p{F\X) dXdF = Yld=i ^d- Both terms in (30) are analytically tractable, with 
the first having the same analytical solution as the one derived in [8]. Further calculations in the the P^ term 
reveal that the optimal setting for q{\id) is also a Gaussian. More specifically, we have: 

:Fv= j q{ud) log — i^^^dud - A (31) 
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where A is a collection of remaining terms and a.d is the mean of (28). (31) is a KL-like quantity and, therefore, 
q{ud) is optimally set to be the quantity appearing in the numerator of the above equation. So: 

q{ud) = e\ ^ ^/^(^)p(urf), 

exactly as in [8]. This is a Gaussian distribution since p{ud) = JV{ud\0^ Kmm)- 
The complete form of the Jensen's lower bound turns out to be: 

D 

Mq.O)= Y,MQ.O)-KL{q\\p) 

^ \^(27r) 2 1/3^2 + ^MMh J 2 2 

O 1 ^ 1 ^ 

- I log \Kt\ - 2 E [Tr i^t-'S,) + Tr (irr V.mD] + 2 ^^''^ 1^^' ^ ^''''^^ ^^^^ 

q=l q=l 

where the last line corresponds to the KL term. Also: 

^0 = Tr{{KNN)qiX)) , ^1 = {KNM)q{X) , ^2 = {KMNKNM)q{X) (33) 

The ^ quantities can be computed analytically as in [8]. 

B Derivatives of the variational bound 

Before giving the expressions for the derivatives of the variational bound (30), it should be reminded that the 
variational parameters iiq and Sq (for all qs) have been reparametrised as Sq = {K^^ + diag{Xq)) and /Xg = 
Kffiq, where the function diag{-) transforms a vector into a square diagonal matrix and vice versa. Given the 
above, the set of the parameters to be optimised is (^/, Ox^ {jlq-, Ag}^^^, X. The gradient w.r.t the inducing 
points X, however, has exactly the same form as for Of and, therefore, is not presented here. Also notice that 
from now on we will often use the term "variational parameters" to refer to the new quantities jlq and Xq. 
Some more notation: 

1. Ag is a scalar, an element of the vector Xq which, in turn, is the main diagonal of the diagonal matrix A^. 

2. Sij = Sq-ij the element of Sq found in the i-th row and j-th column. 

3. Sq = {Sq-ii}fLi, i.e. it is a vector with the diagonal of Sq. 

B.l Derivatives w.r.t the variational parameters 



where: 



J^{q, 6) _ pD l^^o ^^^^TrT 

r/flq 2 r/flq y uf^q 



-Tr 



2 [^/x. 






(35) 
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PD ,^0 ..^^,(^yyT^^^-^ 



Z "^^qiij 



'^^o;i, 



/? 



Tr 






{DK^'m - r'DA-^ - A-'^jYY'^^iA-') 



(36) 



withA = p-'^K_ 



MM 



^2 



B.2 Derivatives w.r.t 6 = {9f, O^) and /3 

Given that the KL term involves only the temporal prior, its gradient w.r.t the parameters Of is zero. Therefore: 

t?7; 'dJ' 



m. 



^0f 



(37) 



with: 






{DK]^'j^ - p-^DA-^ - A-H^YY^-^iA-^) 



^*2 



^e 



/ 



(38) 



The expression above is identical for the derivatives w.r.t the inducing points. For the gradients w.r.t the /3 
term, we have a similar expression: 



D (Tr(if^\^*2) + (iV - M)/?-i - *o) - Tr(yy^) + Tr{ A-'^ J YY^^,) 



;^_ ir 

+I3-^DTi{KmmA-^) + /?-iTr {K]^'j^A-^^JYY'^^iA-^) 



(39) 



In contrast to the above, the term ^y does involve parameters O^, because it involves the variational parame- 
ters that are now reparametrized with Kt, which in turn depends on O^. To demonstrate that, we will forget for a 
moment the reparametrization of Sq and we will express the bound as F{dx, iiq{dx)) (where Hq{dx) = Ktfiq) 
so as to show explicitly the dependency on the variational mean which is now a function of dx- Our calcu- 



lations must now take into account the term ( J-^"' 






that is what we "miss" when we consider 



l^q{Ox) = t^q- 



mx 



l9Ty{ex,IIq{ex)) ^ l}Ty{ex,tlq) f^fyiS] ^ ^l^q{0.c) 

\ '&f^q ) ^Ox 



^Ox 



^Ox 



^Ox 



(40) 



We do the same for Sq and then we can take the resulting equations and replace fiq and Sq with their equals 
so as to take the final expression which only contains jlq and A,: 



^Mex,fiq{ex),Sq{ex)) 

^x 



Tr 



( 



BqKtBq 



IJ-qlJ'q 



{I-BqK.)d^ag(^^^{I-BqK.y 



^Kt 



dOx 



( dP{tlq) V ^Kt 



\ ^fJ-q 



■dOx 



^J-q 



(41) 
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where Bq = A| B~^Aq . and Bq = I -\- Aq KtAq . Note that by using this Bq matrix (which has eigenvalues 
bounded below by one) we have an expression which, when implemented, leads to more numerically stable 
computations, as explained in [9] page 45-46. 

C Predictions 

C.l Predictions given only the test time points 

To approximate the predictive density, we will need to introduce the underlying latent function values F^ G 
R^* ^^ (the noisy-free version of Y^) and the latent variables X^ G R^* ><^. We write the predictive density as 

p{Y,\Y) = f p{Y,,F,,X,\Y,,Y)dF,dX, = f p{Y,\F,)p{F,\X,,Y)p{X,\Y)dF,dX,. (42) 

The term p{F^\X^^Y) is approximated according to 

q{F,\X,) = / n P{^^Md.X,)q{ud)dud = H g(f*,d|X*), (43) 

-^ deD deD 

where q{f^^d\X*) is a Gaussian that can be computed analytically. The termp(X^ |F) in eq. (42) is approximated 
by a Gaussian variational distribution q{X^), 

. Q 

p{X,\Y) ^ / p{X,\X)q{X)dX = (p(X*|X))^(^) = q{X,) = [] ^(x*,,), (44) 

where p{X^^q\X) can be found from the conditional GP prior (see [9]). We can then write 

x*,g = axg + €, (45) 

where a = K^^K^^ and e ~ JV (o^K^^ ~ ^^nk-^Kn^)- Also, Kt = kx{t^t), K^n = ^x(t*,t) and 
i^** = kx{t^t^). Given the above, we know a priori that (44) is a Gaussian and by taking expectations over 
q{X) in the r.h.s. of (45) we find the mean and covariance of q{X^). Substituting for the equivalent forms of 
fiq and Sq from section 2.2 we obtain the final solution 

/ia,^^ = KnPq (46) 

var(x*,g) = h^ - KN{Kt + A^^)"-^kAr*. (47) 

(42) can then be written as: 

p{Y,\Y) = f p{Y,\F,)q{F,\X,)q{X,)dF,dX, = f p{Y,\F,) {q{F,\X,))^^^^^dF, (48) 

Although the expectation appearing in the above integral is not a Gaussian, its moments can be found analyti- 
cally [9, 13], 

E(F,) = 5^^* (49) 

Cov(F*) = B^ (^5 - ^*(^*)T) B^%I- Tr [(i^^\^ - {Kmm + /3^2)"') %] ^, (50) 

where B = p {Kmm ^ ^^^2)'^ "^JY, ^ = (^/(X*,X*)), ^* = (Km.) and ^J = {Km.K.m}- All 
expectations are taken w.r.t. q{X^) and can be calculated analytically, while Km^ denotes the cross-co variance 
matrix between the training inducing inputs X and X^. Finally, since Y^ is just a noisy version of F^, the mean 
and covariance of (48) is just computed as: E(K^) = E(F*) and Cov(F*) = Cov(F*) + I3~'^In^ • 
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C.2 Predictions given the test time points and partially observed outputs 

The expression for the predictive density p{Y^\Yi^^ Y) follows exactly as in section C. 1 but we need to compute 
probabilities for Y^ instead of Y^ and Y is replaced with (F, 1?) in all conditioning sets. Similarly, F is 
replaced with F^. Now q{X^) cannot be found analytically as in section C.l; instead, it is optimised so that 
y? are taken into account. This is done by maximising the variational lower bound on the marginal likelihood: 

p{YP,Y) = J p{YP,Y\X,,X)p{X,,X)dX,dX 

p{Y^ \X)p{YP, YP\X,, X)p{X, , X)dX«dX, 



Notice that here, unlike the main paper, we work with the likelihood after marginalising F, for simplicity. 
Assuming a variational distribution q{X^,, X) and using Jensen's inequality we obtain the lower bound 



/ 



,iX.,X) log K>--|^M>-i';flf;.^M^^.^)dx.dX 



J q{X)\ogp{Y"^\X)dX 

KL[q{X,,X)\\p{X,,X)] 



q{X.,X) 
q{X„X) \ogp{YP, YP\X„ X)dX,dX 



(51) 



This quantity can now be maximized in the same manner as for the bound of the training phase. Unfortunately, 
this means that the variational parameters that are already optimised from the training procedure cannot be used 
here because X and X* are coupled in q{X^,X). A much faster but less accurate method would be to decouple 
the test from the training latent variables by imposing the factorisation ^(X*, X) = g(X)g(X«). Then, equation 
(51) would break into terms containing X, X* or both. The ones containing only X could then be treated as 
constants. 



D Additional results from the experiments 
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(a) 



(b) 



Figure 3: The values of the scales of the ARD kernel after training on the motion capture dataset using the RBF (fig: 
(a)) and the Matern (fig: (b)) kernel to model the dynamics for VGPDS. The scales that have zero value "switch off" the 
corresponding dimension of the latent space. The latent space is, therefore, 3-D for (a) and 4-D for (b). Note that the scales 
were initialized with very similar values. 
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Figure 4: The prediction for two of the test angles for the body (fig: 4(a)) and for the legs part (fig: 4(a)). Continuous line 
is the original test data, dotted line is nearest neighour in scaled space, dashed line is VGPDS (using the RBF kernel for the 
body reconstruction and the Matem for the legs). 







(a) 



(b) 



(c) 



(d) 



Figure 5: Some more examples for the reconstruction achieved for the 'dog' dataset. 40% of the test image's pixels (figures 
(a) and (c)) were presented to the model, which was able to successfully reconstruct them, as can be seen in (b) and (d). 
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